#include <math.h>

#include <iostream>

#include "Interpolation1D.h"
#include "test.h"

typedef double (*Func)(double);

template <class InternalType = double, class ExternalType = float>
void test_subcase(Interpolation1DBase<InternalType, ExternalType> &b,
                  ExternalType xa[], ExternalType ya[], int m, Func f) {
  b.construct(xa, ya);
  for (int i = 0; i < m; i++) {
    ExternalType x = sqrt(i);
    // 1e-6 because float has 6 exact numbers
    ExternalType error = (f(x) == 0)
                          ? fabs(b.evaluate(x) - f(x))
                          : fabs(b.evaluate(x) - f(x)) / f(x);
    CHECK(error < 1e-6);
  }
}

template <class InternalType = double, class ExternalType = float>
void test_case(int m, Func f, int deg) {
  ExternalType *xa = new ExternalType[m];
  ExternalType *ya = new ExternalType[m];
  for (int i = 0; i < m; i++) {
    xa[i] = sqrt(i);
    ya[i] = f(xa[i]);
  }

  // SUBCASE("Polynomial Interpolation") {
  //   // std::cout << "Polynomial Interpolation" << std::endl;
  //   switch (deg) {
  //     case 1: {
  //       PolInt1D<InternalType, ExternalType, 1> p1 =
  //           PolInt1D<InternalType, ExternalType, 1>(m);
  //       test_subcase(p1, xa, ya, m, f);
  //       break;
  //     }
  //     case 2: {
  //       PolInt1D<InternalType, ExternalType, 2> p2 =
  //           PolInt1D<InternalType, ExternalType, 2>(m);
  //       test_subcase(p2, xa, ya, m, f);
  //       break;
  //     }
  //     default: {
  //       PolInt1D<InternalType, ExternalType, 3> p3 =
  //           PolInt1D<InternalType, ExternalType, 3>(m);
  //       test_subcase(p3, xa, ya, m, f);
  //       break;
  //     }
  //   }
  // }

  SUBCASE("Linear Interpolation") {
    // std::cout << "Efficient Linear Interpolation" << std::endl;
    LinInt1D<InternalType, ExternalType> p = LinInt1D<InternalType, ExternalType>(m);
    test_subcase(p, xa, ya, m, f);
  }

  SUBCASE("Quadratic Interpolation") {
    // std::cout << "Quadratic Interpolation" << std::endl;
    QuaInt1D<InternalType, ExternalType> p = QuaInt1D<InternalType, ExternalType>(m);
    test_subcase(p, xa, ya, m, f);
  }

  SUBCASE("Spline Interpolation") {
    // std::cout << "Spline Interpolation" << std::endl;
    SplInt1D<InternalType, ExternalType> s = SplInt1D<InternalType, ExternalType>(m);
    test_subcase(s, xa, ya, m, f);
  }
  delete[] xa, ya;
}

TEST_CASE("Test: Linear Polynomial") {
  int m = 10;
  test_case<double, float>(m, plinear, 1);
}

TEST_CASE("Test: Quadratic Polynomial") {
  int m = 10;
  test_case<double, float>(m, pquadratic, 2);
}

TEST_CASE("Test: Cubic Polynomial") {
  int m = 10;
  test_case<double, float>(m, pcubic, 3);
}

TEST_CASE("Test: Trigonometric Function") {
  int m = 10;
  test_case<double, float>(m, fsin, 3);
}